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In this work we introduce a novel weighted message-passing algorithm based on the cavity method 
to estimate volume-related properties of randonipolytopes, properties which are relevant in various 



Oto estimate voiume-reiatea properties 01 ranaompoiytopes, properties wmcn are relevant m variou: 
research fields ranging from metabolic networks [l[ , to neural networks , to compressed sensing 0] 
■ Unlike the usual approach consisting in approximating the real-valued cavity marginal distributions 

by a few parameters, we propose an algorithm to faithfully represent the entire marginal distribution. 
We explain various alternatives to implement the algorithm and benchmark the theoretical findings 
by showing concrete applications to random polytopes. The results obtained with our approach are 
found to be in very good agreement with the estimates produced by the Hit-and-Run algorithm 
^3 , known to produce uniform sampling. 



I. INTRODUCTION 



There are many problems appearing in various research fields that can be mathematically formalised as having to 
determine the solution set of a collection of linear equalities or inequalities for real- valued unknown variables. More 
^ . precisely, given a rectangular M x TV matrix £ with real entries fx = 1, . . . , M and i = 1, . . . ,N and given a vector 
^ " 7 6 R M , with entries 7 M with fi = 1, . . . , M, one looks for the set V of vectors x <E D C 1^ which are solutions to 
i— —i, e.g. a set of equalities 

(N- N 

^tfa^T", a*=1...-,M (1) 
l^i 

OV 

" or a set of inequalities 

od : N 

E^>7 M , /i=l,...,M. (2) 

(N- *=i 

• • ■ A first example of a problem related to a set of equalities like (JXJ) is Flux Balance Analysis This is a technique 
. \ designed to estimate the reaction rates of a set of chemical reactions in the stationary state. Having in mind future 
applications in this area, we devote some lines to explain this method. Suppose we have a set of N chemical reactions 
' that produce and consume M chemical substances, and let us denote as c M the concentration of the chemical substance 
?3 , fx. Then we can write down the following mass-balance evolution equations 

i— 1 

where xi is the reaction rate of reaction i, 7^ is an exchange flux of chemical substance \x with the environment, and 
(£i i ■ • ■ >£i ) are the stoichiometric coefficients of reaction i. Generally, each reaction rate Xi is a very complicated 
function of the concentrations and kinetic parameters, 1.6. Xi — X 1 (c, kinetic coefficients, etc). Even in cases in which 
all these functional relations are known, the resulting equations are highly non-linear and numerically hard to handle. 
However, if we have good reasons to assume that the system works in the stationary state, that is dc d ^ = 0, then the 
reaction rate vector x can be considered as an unknown vector to a set of equalities ([1]). Besides, due bio-chemical 
constraints, each reaction rate Xi would we such that Xi G [x™ in , x' nax ] so that D = Yii=i[ x i Un > x T ax }- Thus, in this 
context, estimating the volume of solutions corresponds to estimate reaction rates of a set of chemical reactions in 
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the stationary state. Of course, the problem has been transformed from a potentially numerical intractable set of 
coupled non-linear differential equations, to the more tractable yet very time-consuming task of estimating the volume 
of solutions of (|T|) . Flux Balance Analysis goes one step further and simplifies the problem by replacing the whole 
volume V by one of its solutions. This is achieved by imposing an objective function which must be optimized (see 
[l| for details). 

One example for the case of inequalities ((2]) is von Neumann's expanding model for linear economies 0-0] , which, due 
to its interesting applications to metabolic networks i deserves a couple of lines of explanation. In the economic 

context it was originally introduced, this model assumes that an economy is constituted by N companies consuming 
and producing M commodities. Each company i is able to operate linearly with a scale of operation X{. If A = (af) 
and B = (bf) are the input and output matrices of such an economy, as some of the input will be used to produce 
output, the ratio of input to output produced cannot be larger then the a global growth rate p, that is 

P< ^IT 1 T OT $>f-paf)z 4 >0, Vm=1,...,M. (3) 

Z^i=l a i X i i=l 

Thus in von Neumann's model of linear economies, one then wonders what are possible values for the vector of 
operations x for a given growth rate, and how can the optimal growth rate be achieved. 

These two examples are just a glimpse into a myriad of problems, old (e.g. Gardner's optimal capacity problem Q) 
and new ones (e.g. compressed sensing @, [13] ) , coming from diverse research fields which can be mathematically 
formalised as either (TTJ) or (|2|), and which have the same goal: find techniques, either analytic or numerical or a 
combination of both, to estimate the volume of solutions V and statistical properties related to it. 
Thus, rather than focusing in a research field in particular, in the present paper we analyse the problem on general 
terms, and, for sake of clarity, we present our new methodology and test it on simple examples, leaving more complex 
and certainly more exciting applications, particularly in the area of metabolic networks, for future research. 
We have organised this paper as follows: in section [TT] we set up the problem at hand and present it in an unifying way. 
In section UTTl we apply the cavity method to obtain a set of consistency equations for the cavity marginals. We also 
discuss a set of consistency equations for the support of these marginals. As we will see, the set of equations for the 
supports decouples from the actual shape of the marginals, yielding an extremely simple set of equations, which can be 
readily solved. Section HVl contains the first main result, consisting in reweighting the cavity equations seeking for an 
efficient way to solve them. In section [V] we propose two main ways to implement the weighted cavity equations: the 
method of the histograms and the weighted population dynamics in the instance. Besides, for the second method we 
suggest two alternatives: random assignment & variable locking and variable fixing. We have benchmarked all these 
novel ideas with some examples and reported the results in section PVTl Further theoretical research and applications, 
particularly in the area of metabolic networks, are discussed in the conclusions in the last section IvTlI 



II. MODEL DEFINITIONS 



As the analytical treatment to the two problems (p} and © is not that different, we choose to discuss the more 
general problem related to the set of inequalities ([2]). We are particularly interested in heterogeneous systems in which 
the rectangular matrix £ is generally random. This sums up to say quite simply that we focus here in estimating 
volume-related properties of random polytopes. In the so-called H-representation, a polytopc is defined as the set 
of points x = (xi, . . . , xn) £ D = Z>i x ■ ■ ■ x C 1^ encapsulated by M hyperplanes 7 M )}j^Li, with 

^ = (£i, . . . , the normal vector of the /i-hyperplane. Its volume can formally be written as 

From all possible questions related to polytopes in H-representation, we are particularly interested in that of the 
volume V and its projection onto each axis, in order words, the single-site marginal pdf Pi(xi). These two quantities 
can mathematically be written as: 




where we have defined h^xg^) = J2iedp,^i Xi ~ 7^' ®( x ) i s the Heaviside step function, and x^ denotes the vector 
x without component Xi, and D\i = D± x • • • x -D;_i x Z?;+i x • ■ ■ x Djy. Notice that, in a strict mathematical sense, 
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the volume of the polytopc V defined by eq. (j4]) is always strictly 0, as this is the volume of a N — M dimensional 
manifold in N dimensions. This can be formally dealt with as explained in [l3|], but as defined here does not impact 
the results concerning the marginal distributions. 

III. SINGLE-SITE MARGINALS, CAVITY EQUATIONS AND SUPPORTS 



Using the cavity method (see appendix lA]) it is possible to express the single-site marginals Pi(xi) in terms of the 
local variables connected to it. Assuming that the bipartite graph associated to the rectangular matrix £ is locally 
tree-like, we can arrive to the following set of cavity equations 



mf{xi) 



II "^fcOi Vi, v G di, 

7-T / dxg^Q (hnixg^i) + £?Xi) pj^ (xt) , Vi , (1 G di 



(5) 
(6) 



eedfj.\i 



Here the and mft are normalising constants and we have labelled the variable- nodes with Latin indices 
and the factor- nodes with Greek ones Besides, we have used the following standard notation: di denotes the 

set of factor-nodes neighbouring i, dfi denotes the set of variable-nodes neighbouring the factor-node fi, and xg^ is 
the set of dynamical variables on the set variable-nodes dfi. Finally, if A is a set of indices and iei then A\i is the 
set A without index i. 

Once the set of equations ((SJ and ([BJ is solved, the actual single-site marginals are given in terms of the marginals 
{mf{ Xi )}: 



Pi(Xi) 



(7) 



with Vi a normalising constant. It is important to point out that the support of the marginal Pi(xi) is not the 
integration domain for a;, G Di, but rather the result of intersecting the polytope and -D, being the intersection then 
projected onto the ij-axis. Looking at the eqs. ([5]) and ([6]), one notices that it is possible to write down self-consistency 
equations for the support of the marginals in the cavity equations. To do so, let us denote as and Kjy the support 



of the marginals (xi) and P^ v ^(xi), respectively. Then from eq. ([5]) we can write that 



K 



n 



(8) 



To be able to write down the support R^ in terms of k\ v \ we notice from eq. ([6]) that, geometrically, we are 



integrating on a parallelotope S)l ' C Dg^\i defined as the cartesian product of supports {K^}^^, that is 

Xj£edfj,\iKe ■ This parallelotope has vertices, whose set we denote as Vji . As geometrically suggested by eq. 

when Xi is varied the hyperplane + £f Xi = will intersect with all these vertices. Let us denote as 

X ^V^ 1 ^ the set of 2l^'\ 4 l values of Xi at which this intersection occurs. Then the support is the maximal set 
such that 



?(0 



mini ( V 



Notice that once the collection of supports R^ are known, then supports Ki for the single-site marginals Pi{x{) are 
given by 



, max I ( V 



;(') 



(0 

A* 



(9) 



(10) 



The set of equations (JSj> and ((9J is decoupled from the set of the cavity equations ([5]) and ((6j), in the sense that the 
actual shape of the marginals is not needed to obtain information solely about their support. This rather simple 
observation is quite relevant for two reasons. Firstly, to apply the algorithms we devise in Sect. IIVI we need to know 
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the support of the marginals beforehand in order to avoid rejection. This will become clearer in the explanation of the 
method, but in anticipation and with a modest amount of foresight we note that by simply inspecting eq. ([5]): either 

Xi lies within the intersection of all the supports {R$ } fj,edi\v or -at least- one of the m-functions is zero for that value, 
and therefore the value xi should not be allowed. Secondly, as we will see in the derivation below, the self-consistency 
equations for the supports are much simpler than the corresponding equations for the marginals, meaning that they 
can be solved extremely fast. Thus, as the supports give valuable information of allowed values for the dynamical 
variables, these self-consistency equations are important in their own right, as they can be applied, for instance, to 
study the behaviour of metabolic networks under perturbations (e.g. effect of gene knockout in reaction rates) or, as 
we illustrate below in section fVI B[ to determine the critical line for the global growth rate in Von Neumann's model 
described in cq. (|3|). 



A. Expressing the self-consistency equations for the supports in terms of endpoints 

We now move on to implement the set of equations © and © to obtain each support . As the marginals are 
the result of projecting a convex polytope, we expect them to be singly-supported. We then define = [k^l, k^)] 

and similarly R$ = [r„ _,rj*j + ]. From the Heaviside funcion in eq. (|6|) we know that h tl (xg fJi \ i ) + £?Xi > 0. This 
implies that for an auxiliary variable z^ we can write h^{xg^\i)+^ Xi = or Xi = — h^xg^)]. If the reader 

finds the introduction of the variable z M rather obscure, an alternative interpretation will be given in Sect. IIVI Next, 

let us now denote as x 1 ^ ag the 2l 9 ^\ l l values of Xi at each of the vertices of the parallelotope sjx' = x^eSM*^ > 
that is: 

^ £<£dfj.\i 

To look for the maximum and minimum in the set of values {x^,cr a(Ai }, we notice that we do not need to check all 
2I3mVI values, but rather we can write down directly an expression for them. Indeed, let us define first 

(11) 



where z m i n = and z max are the minimum and maximum values the variable z M can take [Hj]. Then the maximum 
and minimum values of 

^ + = max{t« /$ } CTe{ „, +} , <g_ = rrintf^/tf }«={-,+} (12) 



From here we can have finally write 



12) r m m 1 n • r 1 II) T 



(13) 

The set of equations ((TTj) , (|T2|) . and (|T3|) . is a set of self-consistency cavity equations for the endpoints defining the 
supports of the cavity marginals. These cavity equations can be solved in the standard manner by using belief- 
propagation, i.e. fixed-point iteration method, using as a starting support for K\ = Di = [mi_,mi.+] for all 

i = 1, . . . , N and for all fi G di. From cq. (JTUJ) and once the values of have been estimated, we can calculate 

the support Ki of the single-site marginals Pi(xi), viz. 

k it - = maxjr^Ll^Gdi , k it+ = min{r^; + }^ eai 



IV. A WEIGHTED BELIEF PROPAGATION ALGORITHM 



Before presenting our method let us discuss briefly the various alternatives found in the literature of how to estimate 
volume-related properties of polytopes. They can be grouped into two main categories: (i) Monte Carlo simulations 
and (ii) theoretical approaches. 

By Monte Carlo simulations we generally mean a numerical way of directly sampling -hopefully uniformly- the volume 
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of a polytope. The Hit-and-Run algorithm [||, explained in the appendix [B] is a Monte Carlo method that samples 
uniformly the volume of a polytope. This method is reasonably efficient only for small dimensions as its mixing time 
goes as 0(N 3 ). The MinOver + algorithm [3] was originally introduced in neural networks to check whether a solution 
to a problem of the type exists. This algorithm was subsequently applied to sample the volume of polytopes in 
0, UEll- Unlike the Hit-and-Run algorithm, MinOver + can deal with larger systems, yielding samples that become 
more uniform the larger the system is. 

By theoretical approaches what we mean is: one identifies the quantities of interest and under not unreasonable 
assumptions one then tries to write down closed equations for these quantities. This is what we have done here thus 
far (and also in Q) by taking as quantities of interests the single-site marginals, finding the set of closed equations 
(O and © by assuming the Bethe approximation for the problem @. For FBA, that is for problems of the type ([T]) 
this was done in e.g. [l3j . while for the context of compressed sensing similar equations were derived in 0, [l2j|- Of 
course the problem is then how to solve the resulting equations (e.g. set of equations (|5|) and ^ in our case) in an 
efficient way. 

The reader should note that, in the simplest possible analysis, we expect the computational time associated to solving 
the cavity equations to go linear with the system size N (as opposed to ~ N 3 for the Hit-and-Run algorithm), as 
this is precisely how the number of cavity equations grows with N. Leaving this aside there are, however, more 
urgent and pressing matters. Indeed, on general terms it is considered a daunting task trying to solve the exact 
cavity equations by any numerical means when the dynamical variables are continuous, as it is in our case, in 
the cases [3, [H, [l3( or, more generally, for any continued- valued spin models on diluted graphs. To overcome this 
numerical burden, one first approximates the cavity equations somehow and then solves numerically the approximated 
equations. This can be presented in various ways, but the general approach is first noticing that the marginals can be 

parametrised by an infinite number of parameters, that is, (^Xi oif^j an d m$ (xi f3$ where we have denoted 

ol\^ = (a\ V i , j • • • J an( l = f/^i>/^*2> • • ■) ■ For instance, in this framework the Gaussian approximation 
means choosing the parameters of the marginals to be their cumulants and assume that only the first two cumulants 
are different from zero. 

However, for this type of problems consisting in counting the number of solutions of a set of either linear equalities or 
inequalitcs, as it is our case, it is possible to tackle the cavity equations ([5]) and ^ without approximation. To whet 
the appetite for the forthcoming discussion let us first rewrite the eq. © as follows: 

™${xi) = —Try / dzn I dxg^S (h^xg^i) +&Xi -Zfj) (xi) , Vi , p € di (14) 

where we have expressed the Theta functions in terms of a Dirac delta using the integral identity 9(x— a) = J R+ dy8[(x— 
a) —j/]. Note that the new integration variable can be understood as a random variable with uniform density in a 
subset of M + . Indeed, due to the constraints of the problem at hand, the random variable z^ will not generally take 
values on the whole positive line but rather on a compact set of it, and it will be therefore normalisablc. 
Looking at the cavity equation (fT4)) we notice that, written in this way, the Dirac delta suggests to estimate the 
integral expression by using the method of population dynamics (Nb. population dynamics is a smart way to estimate 
the integral using a Monte Carlo method when the integrand is unknown) as in, for instance, the ensemble cavity 
equations for discrete spin models on locally treelike graphs. There are two important differences, though: (i) the 
cavity equation (fLij) is still on the instance; (ii) not all updates Xi <— (z M — (/^(afyAi))/^ suggested by eq. (fT4|) are 

allowed as there may be integration regions such that Xi ^ R/y and therefore the update must be rejected. The first 
difference is unimportant, the second one is not: if we do not find a way to avoid the rejection region, the resulting 
population dynamics performed on the instance will be quite inefficient. 

Fortunately, there is a way to avoid the rejection region completely. To illustrate how, we set aside all notational 
complicacies of the cavity equations and focus on the problem at hand with the following simpler example: 



mix) = — / 

™ J [0,1 



~[dyiPi(yi) 



S \x + ^2 a t yi - 7 , x > 



where Pi{yi) for i = 1, ...,n are some arbitrary pdf with yi £ [0,1]. We note that this integral is zero unless 
x = 7 — X)i=i a iUi — 0- Thus the effective integration region is the one enclosed by the hypercube [0, 1]™ and the 
hyper-plane 7 = X)"=i a iDi- Deciding to do the integral in a specific order we write 

m(x) = — [ dy 1 p 1 (y 1 ) [ dy 2 p 2 {y2) ■■• / dy n p n (y n )S ( x + V] a i y l - 7 ] 
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As we want to evaluate the above integral by Monte Carlo, we need to reweight each density pi on the new region 
R-i{yi,---,yi-i,l)- We write 



Pi{yi\yi,---,yi-i,i) 



PiiVi) 



Wi(yi, ■ ■ .,2/i-i,7) 
Thus the expression for m(x) becomes: 

1 



i(yi,---,Vi-i,i) 



dyzPiiyi) 



K-Hl/i,— )!/i-ii7) 



m(x) 



m 



n 



dyipi(vi\yi, ■ • .,j/i-i,7) 



=1 J T^i(yi,---,yi~i,i) 



~[wi(yx, . . .,3/i_i,7) 



5 I a; + ^ ajj/j - 7 



In this new form the integral m(x) can be evaluated by Monte Carlo without rejection. This is done as follows: 

1. draw yi according to pi(yi\"f), draw y<x according to (2/2I2/1 » t)> e ^ c i an d finally draw y n according to 

Pn(y n \yi,---,y n -i,i) 

2. Assign x <- 7 - 5Z"=i a i2/i witn weight fl = JX"=i Wi(yi, . . . , j/i-i, 7) 

This process can be then repeated A/J, times to obtain a collection of pairs {(x a , fJa)}a=i,....AT p ; from which we can 

reconstruct m(x) as m(x) sa ^ So^i ^a^^— x a ) with m = JZ^i ^a- This reweighting technique avoids the rejection 
region allowing us to estimate the involved integrals much more efficiently. 
Thus, if wc apply this reweighting technique to cavity equations (|14j) . we obtain 



mf{ Xi ) = — 



n 



dxt ] Pf M) {xi 3 \x^ , . . . , xi>._ 1 , 7 M ) 



(*/!,.••.**,•_! .7") 



x 8 (M^/A'i) + tfxi - Zfi) 



(15) 



where we have used the notation fc^ = \ i| and S ctyi \ i for j = 1, . . . , fc^, and wc have defined 



w,- (a;^,. • .,a*, ,,7^) 



ujf\x £l ,. . .,x^_ 1 ,7 M ) 



dx^f{x t] ) 



Kf (^ 1 ,... > ^._ 1 ,7") 



Although wc could have chosen any order of integration to write down eq. (|15p . it is convenient, as we will explain 
below, to integrate the variable the first one [l^- We will henceforth generally call the set of equations (fTS")) 
and ([5]) as i/ie weighted cavity equations and the algorithm to solve them weighted belief-propagation or weighted 
message-passing algorithm. 



V. IMPLEMENTING THE WEIGHTED BELIEF-PROPAGATION ALGORITHM 

We move on to discuss the actual implementation of the algorithm to numerically solve the set of equations ([5]) and 
(fi3|) . As we will see shortly, although we have found a nice method to avoid rejection, further complicacies appear 
in the horizon which may put into question the discussion we have had thus far. Fortunately, we describe a way to 
overcome them. 

We present here two alternative implementations: the method of histograms and the weighted population dynamics in 
the instance. 
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A. Method of histograms 

Ideally, to solve numerically the weighted cavity equations we would firstly represent the marginals P^ (xi) 
and m^'(xi) with populations of pairs { ( , x^ ) \ and I { v \^l, x $u) \ , respectively, so that 

pM(xi) « (l/P^E^itfg^Xi-xg) and m«(.xO « {1/mf) £^ v$ a 8 (x t - x$ a ) , where and mf 

are normalisation factors. There is a drawback though: while it is possible -and desirable- to use the population of the 
P's to evaluate eq. (|15l) updating the population for the m's, it is difficult to see how to use directly this representation 
to update the population for the m's according to eq. ([5]). Of course, one could always construct histograms for the set 
of functions {m£ i (a;t)}>e8»/v to calculate a histogram for Pj> v '(xi) according to eq. ([5]), and then draw a population 

for i 3 / (xi ) from its histogram. But if we need to construct histograms at some point in the algorithm, it is actually 
a waste of time to go back and forth between a representation of histograms and a representation of populations. 
Thus, due to equation (|5|), we are unfortunately obliged to implement a weighted message-passing algorithm where 
the functions P ^ (xj) and m^(xi) are quite simply represented by histograms. We call this implementation the 
method of histograms. 

Discretisation of the unweighted cavity equations and their Fourier transform was the method used in [l2l [l3j in the 
context of FBA and compressive sensing. Let us see how it is possible to avoid discretisation altogether. 



B. Weighted population dynamics in the instance 



Fortunately, it is possible to overcome the previous problem of having to go through the m-functions via equation 
([5]). To illustrate the method we again ignore all tediousness regarding notation and analyse the following simpler 
example that captures the complicacies arising from combining equations (|5|) and (|15[) . namely, the appearance of 
multiple Dirac deltas: 



r(x) = — 



[o,iv 



~[dyipi(yi) 



[x + V a iyi 7 / 

V *=i / J ^' 1 



~[dziipi(zi) 



six + ^bi 



with x > and where Pi{yi), i = 1, ■ • ■ , n and tpi(zi), i = 1, . . . , m are some arbitrary pdf with y.;, zi € [0, 1]. As we 
have already discussed the problem of rejection, we first simply reweight this equation accordingly: 



rix) 



n 



dyipi(yi\yi, ■ ■ .,2/i_i,7) 



n 



=i«r(»i !«-i.t) 

dZiipi(zi\zi,...,Zi-i,5) 



i=l 



lwf\z 1 ,...,z i . 1 ,S) 



+ ~ T 

i=i j 

m \ 



(16) 



As there are two Dirac deltas in the expression (|16|) . it is natural to ask how resolve the issue of assignment for the 
variable x so that we can estimate the integrals by Monte Carlo. We describe two options: random assignment & 
variable locking and variable fixing. 



1. Random Assignment & Variable Locking 

In this case we use one of the Dirac deltas to assign a value to x, locking such value in the other Dirac deltas. 
Suppose, for instance, that in our example we use the first Dirac delta to assign a value to x. The first steps of the 
algorithm would be almost as before: 

1. draw yi according to pi(yi\j), draw y% according to P2{y2\Di,l), etc, and finally draw y n according to 

p n (y n \yi,---,y n -i,i) 

2. Assign x <- x* = 7 - Yh=i a iyi 

3. Calculate the weight fl = n™=i w \ P \yii ■ ■ • j 2/i-ii7) 

In the second Dirac delta, the value of x has been locked to a;*. The integration regions will depend on x* and we 
denote this by writing (zi , . . . , Zi-i, 7, x*). The multiple integral over z is then done as follows: 
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1. draw Z\ according to f/'itzil^, draw Z2 according to ifai&^z-i, 5, x*), etc, and finally draw z m _i according to 

1pm-l( z m-l\ z l, ■ • ■ 7 z m-2, S, X*) . 

2. The Dirac delta locks the value z m <— (<5 — x* — biZi)/b m allowing us to integrate the last integral over 

3. Calculate the weight T = ip m (z m \zi, . . . , z m -i, S, x*) J]™ 1 wf^ (zi, . . . , 6, x*) = 
ip m (z m ) YlTJi w i ( z u ■■■,z i - 1 ,6,x*) 

4. Assign x* an overall weight A = Or 

As before this process can be repeated Af p times to obtain a collection of pairs {{x^, A a )} a=1 t ...jj > from which we 
can reconstruct r(x) as r(x) sa i So^i A a S(x — x* a ) with r = Y^a=i ^a- 

2. Variable fixing 

In this case we fix the value of x from the beginning, then evaluate each integral according to the second part of 
the algorithm before. For our example at hand the steps are: 

1. Fix a value of x to x* within its support for r(x). 

2. For the integrals involving function p: 

(a) draw yi according to pi(yi|7,x*), draw y2 according to P2 (z/2 ^ li x *)> etc, and finally draw y n -i according 
to p„_i(y„_i|yi,...,y„_2,7,a;*). 

(b) The Dirac delta locks the value y n <— (S — x — a i z i)/ a n allowing us to integrate the last integral over 
z n . 

(c) Calculate the weight T p = ^ n (y n \y 1 ,...,y n _ 1 ,j,x*)l\ r - =1 wl P \yi,...,y i ^i,~f,x*) = 
Pn(y n ) Wt=i w t P) (yii ■■■,Vi-Ul,x*) 

3. Repeat same process for function if) and calculate the weight 

4. Assign a weight A = TpT^ to the fixed value x* 

This process can be repeated M p times, scanning the support of r(x) to obtain a collection of pairs {(x a , A Q )} Q=l! ...jj , 
from which we can reconstruct r(x) as r(x) « ^ A Q £(;r — x a ) with r = J2^=i A<*- 

This implementation has two clear advantages with respect the previous one. The first one is that since we are not 
using the Dirac delta to randomly assign a value of x 7 this can be used to perform one less integral. For this reason, 
we have found more convenient to arrange the order of integration in cq. (|I5|) . so that the integral over z^ is the 
last one to be estimated by the algorithm, so that the contributing weight is the unity. The second advantage lies 
on the fact that at a fix value of x what we should obtain is an averaged value of the weight. The weight is itself 
a random variable whose variance, even in simple cases, can be rather large. Thus it is appropriate to refine the 
estimates of the weights by replacing them by averaged values. To do this we simply calculate T estimates for each 
weight {T pt } t =i t and {r^} t=1 ^, calculate the average weights T p — ^Y^t=i^p,t an d T^, = ^ Ylt=i ^iM an< ^ 
assign instead the weight T p T^, to the fixed value x* [20| ■ 

3. Implementation to our problem 

Considering all the points discussed previously, while hoping that the tedious but necessary notation we are about 
to use does not distract the reader, the core of the algorithm for the weighted population dynamics to solve the cavity 
equations ([5]) and (fT5)) is as follows. We assume that all supports for each marginal P^\xi) are known and 

we represent each marginal by a population of pairs I \ w[^,x\^) \ . Then for the random assignment & 

IV ' ' / J a=l,...,J\f p 

variable locking the essential steps are the following: 

1. Choose a variable-node i, a factor- node v £ di, and a value Xi £ k[ v ^ 

2. Chose a factor- node /jq G di/v 
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(a) Draw Xi t with probability P^ 10 ^^ | 7 Mo ), draw X£ 2 with probability P^"\x£ 2 \xi 1 ,j tJ, °), etc, draw xe (i) 
with probability P^ ^ {xi 1 \xi 1 , ...,xi fi , , 7 M0 ), and draw a uniform random variable z p in the segment 

<°W\i>7 w ) 



(b) Assign x t <- [z Mo - M x %o\i)]/£f 

(c) Calculate 0$ = wMfo^, 7 M ) • • • 7") 
3. For all /i e di/{v\J fi ) 



(a) Draw with probability P^ (x^ 7' 1 ), draw xi 2 with probability P^\xi 2 \xi 1 ,Xi,j fl ), etc, draw (i) 
with probability P/7 (x^ |x fl , . . . , xi ... , x 4 , 7 M ) 

(b) Calculate = ]\% Jf (x £l x tj _ i: x h 7") 

4. Calculate t\"^ = Yi^ai/^^^ ano - replace a P a h of the population of the marginal pj; u \xi) by the new pair 

(ifU) 

For the second implementation, that is variable fixing, we enumerate the following basic steps: 

1. Choose a variable- node i, a factor-node v £ <9i, and a value x^ S -f^"^ 

2. For all [i^dijv 

(a) For t = 1,...,T 

i. Draw x^ with probability P^ (x^ \xi, 7^), draw X£ 2 with probability P^ ( x ^2 l^i 1 x ii 7 M )i e t c i draw 



x<> with probability P/ ^ (x^ |x fl , . . . , xt ... , a;*, 7 

ii. Calculate w t = rij=i w } (^i > • • • > ^-1 > ^t> 7**) 
(b) Set n<? = i ELi «* 

3. Calculate = IlueOi/i/ ®ff anc ^ replace a pair of the population of the marginal pj; u \xi) by the new pair 



f 1 



(ifU) 



It is important to notice here that in the last step the replacement must be done according to how the value of Xi 
is selected on the first place, otherwise we could generate a biased sampling of the support -K.- , that is, if in the 
first step Xi is chosen uniformly randomly, or fixed then the replacement must be done in the same manner. Besides, 
instead of replacing other alternatives would be, for instance, to increase the populations in those regions where the 
sampling is rather poor, or to use a non-uniform sampling taking more points wherever they are most needed. 

VI. NUMERICAL TESTS 

To illustrate the previously discussed results we apply them to some simple examples. 

A. A toy example regarding supports of marginals 

Consider the polyhedron described by the following set of inequalities £ e cc > 7 with 

* = \y \ . 7= (]z e ) a?) 
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Figure 1: Bipartite graphs associated to the polyhedron defined in our example (|17[) . The left one it is exactly a tree corre- 
sponding for e = 0. The right one corresponds to the generic loopy case of e ^ 0. 



For e = the associated graph to £ e is exactly a tree, as shown in figure Q] so that iteration of the cavity equations for 
the support converge at one step. On the other side, if e 7^ the graph is loopy and, depending on the value of e the 
estimates from the cavity equations can be very bad indeed. Yet for small e one obtains very precise results. In figure 
[5] we plotted the results of iterating the cavity equations (fTTj) , ([T12]) . and (fT5|) for the supports. Due to symmetry of 
the problem, the support and marginals is the same for the three variables and we parametrise the support as [a, b]. 
We have taken a value of e = 0.05. 

It is important to notice that as the dynamical variables are continous we have at our disposal an ample set of 
transformations we can perform on the system. This can be used to transform a very loopy network into a more 
tree-like one. For instance, in this very simple example we notice that we can write £ c = rj e T e with 
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Then, the transformation y = T e x takes the system into an exact tree. Of course, under which general conditions 
it is possible to perform these type of transformation is not clear to us just yet, nor it is clear which are the most 
appropriate transformations, whose inverse transformation can be performed efficiently. 
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Figure 2: Results of iterating the cavity equations (|12|l . and (|13|l for the supports. Due to symmetry of the problem, the 

support and marginals is the same for the three variables and we parametrise the support as [a, b]. We have taken a value of 
e — 0.05. Left (right) figure corresponds to the evolution of the cavity equations for the lower endpoint a (b). The inset reports 
the absolute error between the value of the endpoint at iteration t and the exact value. 
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B. The critical line p(n) in Von Neumann's model 

As a second example we revisit the Von Neumann model ([3]). We are interested in calculating the value of the 
critical line p c (n) (with n = N/M) above which it is not possible to find more solutions to the set of inequalities [f|. 
To find this critical line using the self-consistency equations for the supports, we assume that the critical line occurs 
when at least one support becomes zero. Thus for a fixed value of the ratio n = N/M we increase the value of the 
global growth rate p until at least one of the supports becomes zero. 

Numerical results for the critical line are reported in the left panel of figure [3] and compared with the MinOver + 
algorithm. Here we have generated Poissonian graphs both for factor and variables nodes. Factor nodes have an 
in-degree and out-degree average of 3. Thus the variable nodes have the in- and out-average degree equal 3/n. The 
sizes of the graphs have been generated as follows: for n <= 1, M = 500 and N is varied while for n > 1 N = 500 
and M is varied. The results are average over 50 graphs. Although the cavity equations for the supports can cope 
with much larger graphs we keep their sizes small so as to compare with the MinOver + algorithm. 



p c (n) 




0.5 1 1.5 2 2.5 3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

n p/p c 



Figure 3: Left: results for the critical line p c (n) from Minover + (red) and belief propagation for the supports (green). Details 
are explained in the text. Right: plot of the minimum support size (black-dotted lines) and average support size (red-solid 
lines) versus p/p c for the ratios n = 0.5, 2. 

In the right panel of figure [3] we have plotted two possibles "order parameters" that could be used to locate the 
critical line p c . These are the minimum support size A m!Il and the average support size A avg . It is interesting to 
notice that for n < 1 both parameters become zero at the critical line (as given by the MinOver + algorithm). This 
implies that at the critical line the volume of the polytope collapses to a point. Remarkably, for n > 1 we have that 
A m j„ = and A avg > at the critical line, indicating that the polytope has collapsed in one or more dimensions, but 
it is still dimensionfull in a subspace. 

In figure 0] we also report the average running time t corresponding to figure |3] for both the belief propagation of the 
supports and the MinOver + algorithm. As it can be clearly observed, belief propagation is faster than MinOvcr + by 
an order of magnitude. Besides, the average running time i grows with n = N/M for n < 1, while remains constant 
for n > 1. This is not unsurprising considering how the networks are generated (see description in the text above). 



C. A toy example for Weighted Population Dynamics in the Instance 

In this section we illustrate random assignment & variable locking and variable fixing in a particular case of the 
example (| 16[) discussed before. As we want to illustrate the method to estimate the integral rather than its use to find 
the marginals, we consider here the the marginals to be known Gaussian distribution gi(y; rrii, of) = gi{y) normalised 
on the interval [0, 1] and with mean value rrii and variance of. 

r(x) = -I 2 (x), I(x)= f 



~[dyigi(yi) 



_i=l 



3 

i=i 
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0.5 1 1.5 2 2.5 3 



n 

Figure 4: Average running time i versus n for the belief propagation for the supports and the MinOver + . As we can see the 
belief propagation is an order of magnitude faster. 



where we also assume that x G [0, 1]. Although this a pedagogical example, it may be worth it, for the sake of clarity 
-but with the risk of being a bit repetitive- to discuss it in some detail. We describe here how the method of random 
assignment & variable locking looks like on a pen-and-paper calculation. Firstly, for one of the multiple integrals 
I(x), we note that due to the assignment x <— 1 — ^ i=1 2/i, and the fact that x, y%, 2/2, 2/3 G [0, 1], we have that from 
the integration region [0,1] 3 in the (2/1, 2/2, 2/3)~ s Pace, the region that actually contributes to the integral is the one 
below the plane y 3 = 1 — y± — 2/2- Choosing an order of integration we write 

pl ri—yi /-l— 3/1—3/2 I 3 

I(x)= dyxg 1 {y 1 ) I dy 2 g 2 (y2\yi) dy 3 g 3 (y 3 \yi, 2/2)^2(2/1)^3(2/1, 2/2)6 ix + Y\i - 1 

Jo Jo Jo y i=1 

where we have also reweighted the pdfs g 2 and g 3 so that is it normalised in the new integration intervals [0, 1 — 2/1] 
and [0, 1 — yi — 2/2] respectively, with weights 

1/1 

w 2 (yi) = / dy 2 g(y2) = Q(yi\m2,ai), w 3 (yi, 2/2) = Q (2/1 +y2\m 3 ,ai) 
Jo 

where we have defined the function 



Q(x\a, b) = 




and with (72 (j/2 |yi ) = w 2 {yi) ' 52 (2/3 1 2/2) = mffaf 3 ^) ■ Then, to estimate the integral /(a;) and assign a value of a; we 
must implement the following steps: 

1. Draw a Gaussian random variable y\ with mean mi and variance in the interval [0,1], draw a Gaussian 
random variable 2/2 with mean m,2 and variance erf in the interval [0, 1 — y±], and draw Gaussian random variable 
2/3 with mean m 3 and variance erf in the interval [0, 1 — 2/1 — 2/2], 

2. Assign x the value x = 1 — Si=i 2/i 

3. Calculate the weight (2/1)^3(2/1, 2/2) 

Now for the second /(x) in our example, the value of x is now locked at x* . After restricting the integral to the region 
with non-zero value and using the Dirac delta to integrate over z 3 we can write 

pl — x* pX—Zi—x* 

I(x*) = / dz 1 g 1 (z 1 ) / dz2g 2 {z 2 \z 1} x*)g 3 [z 3 (z 1 ,Z2,x*)}w 1 (x*)w2(yi,x*) 
Jo Jo 
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with Zs{zi, Z2,x*) — 1 — x* — z\ — Z2, and with the weights with the same definition as before. From here we see that 
to estimate I{x*) by Monte Carlo we can do the following: 

1. Draw a Gaussian random variable z\ with mean m\ and variance a\ on the interval [0, 1 — x*\ and draw a 
Gaussian random variable z-i with mean m-i and variance a\ on the interval [0, 1 — x* ~ z{\ 

2. Caculate the weight g^z^zi, z 2l x*)]wi(x*)w2{yi 1 x*). 

From here the total weight is u — W2 (2/1)103(2/1, 2/2)273 [^3(21, Z2, x*)]wi(x*)ui2(yi, x*) at x = x* . Repeating the process 
M times we obtain a population of pairs {(s*,w Q )}^ =1 from which to construct r(x). 

To implement the numerics we have chosen mi = 0, 7712 = 1/2, and 772.3 = 1 while for the standard deviations we have 
taken a\ = VOA, a 2 = \/0?2, and 0-3 = \/~0A and we reported the results in figure [5j with a more detailed explanation 
in its caption. 




Figure 5: In this figure we report the results of the weighted population dynamics in the instance applied to the toy example 
reported in IVI CI We have taken 3 Gaussian distributions (left panel) with means 0, 1/2, and 1 and variances a 2 — 0.4,0.2, 
and 0.4. The centre panel corresponds to estimating r(x) by WPDI using the method of random assignment & variable locking 
using a population size of AT = 7 ■ 10 2 , 7 • 10 3 , and 7 • 10 5 (top to bottom). The right panel corresponds to WPDI using the 
variable fixing method for a population of TV = 8 • 10 2 and using averaged weighted using sizes T = 15, 60, 300 (top to bottom). 
The orange solid line in the centre and right panels corresponds to the exact result for r(x) 



D. Random graph and Comparison with Hit-and-Run algorithm 

Finally our last example corresponds to a random network of N = 25 variable nodes and M = 10 factor nodes. 
Although small, this size can be found in actual applications in metabolic networks as, for instance, human red blood 
cells [13, EH- Here we have used WPDI with variable fixing with a population size of M = 10 3 . For the averaged 
weights we use a time window of T = 10 2 weights for the cavity marginals and T = 10 6 for the real marginals. 
Besides, for the final marginals a population is not used due to the simple fact that then we need to plot them, so 
we construct directly the histogram by fixing the value of x at the midpoint at each bin. Results for each marginal 
Pi(xi) for i = 1, . . . , 25 are presented in figure [5] and we have compared our findings with results obtained with the 
Hit-and-Run algorithm (for implementation see appendix [B]) . As one can see the agreement is excellent considering 
that the network is small and therefore loopy and our equations rely on the Bethe approximation. 

VII. CONCLUSIONS 

As we have discussed, many problems arising from diverse research topics may be mathematically recast as properties 
related to the volume of polytopes. Methods to estimate these properties can be roughly divided into mathematical 
and Monte Carlo approaches. In both cases, the approach used may not be efficient enough to treat polytopes of 
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Figure 6: Results for the marginals Pi{xi) for i = 1, ... ,25 (left to right and top to bottom) of the WPDI algorithm with 
variable fixing (blue) compared with estimates for the marginals using the hit-and-run algorithm (red) 



practical relevance (e.g. large metabolic networks). In this paper we have shown that for diluted random polytopes, a 
novel weighted belief propagation algorithm can be used to calculate efficient single-site marginals. We have discussed 
several options of implementing this algorithm, devoting more lines to explain in detail the WPDI with its two versions: 
random assignment & variable locking and variable fixing. Importantly, we also pointed out that it is possible to write 
down self-consistency equations for the supports of the marginals which, apart of its use in implementing WPDI 
efficiently, can be employed to obtain relevant information about the system with minor effort (e.g. to estimate the 
ranges of reaction rates allowed by stoichiometry; to evaluate the impact of a genetic mutation). Examples are used 
throughout this paper to illustrate, to support, and to benchmark our novel theoretical findings and as so, we have 
kept most of them fairly simple. The attentive reader would have noticed that we have left some computational issues 
unattended. For instance we have not discussed in detail the fact that eq. ( Hp is quadratic in the number of neigbours 
of i and how this dependence can be linearised by adapting the nice trick in [13[ . Our main goal here is to provoke the 
reader with the novel ideas presented here rather than discussing ways to improve the algorithmic implementation. 
These issues are certainly important, but we believe they belong in another paper. 

The work presented in this paper has several applications and extensions. First of all we aim to apply the weighted 
belief propagation to real metabolic networks, possibly starting from the human red blood cell, which is of size 
comparable to the example presented here, and extending the study to other reference metabolic networks (e.g. E. 
coli). Since we are able to tackle both FBA and the Von Neumann frameworks, we may use our tools for instance to 
compare these two approaches and to relate our results to the ones already published by using these two methods. 
Secondly, we plan to apply the weighted belief propagation to solve a metabolic network model where the direction of 
the reactions is not fixed and has to be determined by minimization of the Gibbs free energy. Finally, we plan to see 
how to use the self-consistency equations of the supports to study for instance the behaviour of metabolic networks 
under perturbations. 
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Appendix A: Derivation of single-site marginals using the cavity method 

The first step in deriving the cavity equations (O,© for problem (|2J) is to start by calculating the joint marginal 
Pn(%du,) °f the set of variables xg^ connected to the factor node \x. By definition this is given by 



with f^(xg^) = Q[h^(xg^)]. Here, Pjf^(x ail ) is the joint marginal for the variables xg^ in the absence of factor node 
fi. If the graph is a tree or it is locally tree-like, the set of variables xg^ in the absence of fx are mostly uncorrelated, 

and we can confidently write P^\xg^) = Yle^g^ P^\xe)- 

Let us move on to find an expression for the single site marginals Pj(xj): 

Pi(xi) = y ( dx\i JJ/^(x a/1 ) 

= y j dx\i JJ f^Xg^) Y[ U(Xdn) 



V 1 



(d/j,3i)\i 



Here the notation xr afl3 ^\i stands for the collection of variables nodes which share a factor node with i with the 

exception of the node i itself. The important point to note here is that J dxg^i Il^di f^{xa^) oc P^m^ 1 {x(dn3i)\i)- 
As variable node i is absent in this marginalisation and as we are considering locally tree-like graphs we can confidently 
write P^ 9 i {x(d ll3i )\i) = ILeO* p l { x dn\i)- All in all we have the following expression: 



Recalling that the joint marginals Ppf' themselves factorise we finally write 



fj.£di j£dfi\i 

Close equations for the cavity marg inals P^\xj) can readily be written down by simply rcmovoing one of the factor 
nodes in the neighbourhood of i, that is: 

p^Xxi) = i n J dx afl \iU(x aii ) n p^ix,). 

Notice that the functions rn^\xi) we introduced in the main text are the combinations of terms in from of the first 
product in the predecing equation: 



j£dfj.\i 



Appendix B: The Hit-and-Run algorithm 



As we mentioned previously, the Hit-and-Run algorithm is a numerical method to sample directly and uniformly 
the volume of a polytopc. This useful feature makes it particularly attractive for problems such as FBA [l6|], but has 
the unfortunate drawback of being rather slow. Nevertheless, we find it useful in our case as a benchmark for our 
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message-passing algorithm. 

We explain the Hit-and-Run algorithm by following the nice work in (l7| , simply adapting the notation. Consider the 
following set of inequalities 

(efx>r, /i = i,...,M 

with x £ M. N and = 1, V/x = 1, . . . , M the unit vectors normal to the M hyper-planes encapsulating the polytope. 
Leaving aside some mathematical technicalities and simply assuming the polytope is properly defined, the algorithm 
samples the volume by starting at a given point X inside the polytope and moving along a random direction v. To 
know how far to move one needs to calculate all the intcrcctions between all hypcrplancs and the straight line passing 
through X in the v direction. Then the closest hyperplanes would be the furthest region within the polytope that 
can be reached. This procedure can be iterated randomly as follows: 

Step Find a point X^ ' interior to the polytope. Set n = 0. 

Step 1 Generate a direction vector v^ n ' from a uniform distribution. 

Step 2 Determine the intersections 



_ ( g M)T X (n)_ 7 M 

A+ = min {A M |A^ > 0} , A" = max - {A M |A M < 0}, 

where X* 1 allow to determine the closest hyperplanes to X^ n > . 
Step 3 Generate a uniform random number u £ [0,1] and set 

X (n+l) = x (n) + ( A - + U ( A + _ v (n) 

Step 4 Set + 1 and go to Step 1 or stop if satisfied with the sampling. 

It was proven 0, [l7| that the statistics of the sequence {X^}^ =0 converges to the uniform distribution on the 
polytope. 



Acknowledgments 

IPC would like to thank Lenka Zdeborova for pointing out some work done in the area of compressed sensing. FFC 
would like to thank the hospitality of the department of Mathematics at King's College London. FAM aknowlcdgcs 
European Union Grants PIRG-GA-2010-277166 and PIRG-GA-2010-268342. FFC aknowledges Spanish Government 
grant FIS2009-09508. We would also like to thank the referees from JSTAT to help us make this a better well-rounded 
paper. 



[1] K. J. Kauffman, P. Prakash, and J. S. Edwards, Current opinion in biotechnology 14, 491 (2003). 
[2] B. E. Gardner and Derrida, J. Phys. A 21, 271 (1988). 

[3] F. Krzakala, M. Mezard, F. Sausset, Y. Sun, and L. Zdeborova, Preprint arXiv:1109.4424 (2011). 
[4] R. L. Smith, Operations Research pp. 1296-1308 (1984). 
[5] J. V. Neumann, Ergebn. eines Math. Kolloq. 8, 73 (1945). 
[6] A. D. Martino and M. Marsili, JSTAT 2005, L09003 (2005). 

[7] A. D. Martino, C. Martelli, R. Monasson, and I. Perez Castillo, JSTAT 2007, P05012 (2007). 

[8] A. D. Martino, C. Martelli, and F. A. Massucci, Europhysics Letters 85, 38007 (2009). 

[9] C. Martelli, A. D. Martino, E. Marinari, M. Marsili, and I. Perez Castillo, PNAS 106, 2607 (2009). 
[10] A. D. Martino and E. Marinari, in Journal of Physics: Conference Series (2010), vol. 233, p. 012019. 
[11] A. D. Martino, D. Granata, E. Marinari, C. Martelli, and V. V. Kerrebroeck, J. Biomed. Biotech. 2010 (2010). 
[12] D. Baron, S. Sarvotham, and R. G. Baraniuk, Signal Processing, IEEE Transactions on 58, 269 (2010). 
[13] A. Braunstein, R. Mulet, and A. Pagnani, BMC bioinformatics 9, 240 (2008). 
[14] W. Krauth and M. Mezard, J. Phys. A 20, L745 (1987). 

[15] D. D. Martino, M. Figliuzzi, A. D. Martino, and E. Marinari, Preprint arXiv:1107.2330 (2011). 



17 



[16] K. Almaas, B. Kovacs, T. Vicsek, Z. M. Oltvai, and A.-L. Barabasi, Nature 427, 839 (2004). 

[17] H. Berbee, C. Boender, A. R. Ran, C. Scheffer, R. Smith, and J. Telgen, Mathematical Programming 37, 184 (1987). 
[18] Since the polytope is at most restricted by D, then z M has indeed a finite maximum value. For practical purposes one 

simply puts a cutoff maximum value. 
[19] Note that while is the variable to formally be integrated the first one, within our Monte Carlo method, z M is the last 

variable to be estimated. 

[20] Clearly, T does not have to be fixed, but chosen so as to achieve certain accuracy, nor does it have to be the same for all 
weights. 



